---
title: "Lab 4 Code"
author: "Fabio Setti"
date: "Spring-2026"
output: html_document
---

## Today's Packages and Data 🤗



```{r}
library(tidyverse)
library(car)
theme_set(theme_classic(base_size = 16, 
          accent = "#3492eb", 
          base_family = 'serif'))
```


```{r}
dat <- rio::import("https://fabio-setti.netlify.app/data/World_happiness_2024.csv")
```


## Partial regression coefficients


```{r}
reg_full <- lm(Happiness_score ~ Corruption + Freedom, dat)
coef(reg_full)
```



```{r}
reg_free <- lm(Happiness_score ~ Freedom, dat)
coef(reg_free)

reg_corr <- lm(Happiness_score ~ Corruption, dat)
coef(reg_corr)
```

## Partial regression coefficients "By Hand"


```{r}
resid_happy_Free <- resid(lm(Happiness_score ~ Freedom, 
               dat))

resid_corr <- resid(lm(Corruption ~ Freedom, dat))

# the 0 +  the intercept, which is exactly 0 anyway in this case 
coef(lm(resid_happy_Free ~ 0 + resid_corr))
```



```{r}
resid_happy_corr <- resid(lm(Happiness_score ~ Corruption, 
               dat))

resid_Free <- resid(lm(Freedom ~ Corruption, dat))

# the 0 + removes the intercept, which is exactly 0 anyway in this case 
coef(lm(resid_happy_corr ~ 0 + resid_Free))
```


## Added variable plots 


```{r}
ggplot() +
  geom_point(aes(y = resid_happy_Free, 
                 x = resid_corr))
```


```{r}
ggplot() +
  geom_point(aes(y = resid_happy_corr, 
                 x = resid_Free))
```


## The Quick way: `car::avPlots()`   


```{r}
car::avPlots(reg_full)
```


## Another Look at the 3D plane


```{r}
library(FabioFun)

nice_3D_plot(y = dat$Happiness_score,
             x1 = dat$Corruption,
             x2 = dat$Freedom,
             axis_names = c("Happiness", 
                            "Corruption", 
                            "Freedom"),
             reg_plane = TRUE)
```


## "By Hand" Example of Bootstrapping

```{r}
# empty element to save R^2 to
r_squared <- c()

reg_vars <- dat[, c("Happiness_score",
                        "Corruption",
                        "Freedom")]

set.seed(34677)

for(i in 1:2000){

# sample from data 
sample <- sample(1:nrow(reg_vars), 
                 replace = TRUE)  

dat_boot <- reg_vars[sample,]

# Run regression
reg_boot <- lm(Happiness_score ~ Corruption + Freedom, 
               dat_boot)  

# Save R^2
r_squared[i] <- summary(reg_boot)$r.squared
}
```



```{r}
summary(reg_full)$r.squared
```


```{r}
hist(r_squared)
```


```{r}
quantile(r_squared, c(.025, .5, .975))
```


## The Quick way: `Boot()` function from `car`  


```{r}
set.seed(875)
boot_model <- car::Boot(reg_full, R = 2000)

confint(boot_model, type = "perc")
```


```{r}
confint(reg_full)
```

## Visualizing Bootstrap Results


```{r}
hist(boot_model, 
     parm = "Corruption")
```


```{r}

hist(boot_model, 
     parm = "Freedom")
```